#include "util.hpp"
#include "qual.hpp"

// TODO: in all quality functions the read_length and original_bin_count should be queried from the loaded matrix.

void rebin_quals(const array_type5 &calls, unsigned int read_length, unsigned int orig_bin_cnt, unsigned int new_bin_cnt, boost::multi_array<int, 2> &new_bins)
{
    boost::multi_array<double, 2> accum(boost::extents[read_length][orig_bin_cnt]);
    new_bins.resize(boost::extents[read_length][orig_bin_cnt]);

    unsigned int i,j,ii,jj,kk;

    for(i = 0; i < read_length; i++)
        for(j = 0; j < orig_bin_cnt; j++)
            for(ii = 0; ii < 4; ii++)
                for(jj = 0; jj < 4; jj++)
                    for(kk = 0; kk < 4; kk++)
                        accum[i][j] += calls[i][j][ii][jj][kk];
    
    for(i = 0; i < read_length; i++)
    {
        double sum = 0;
        for(j = 0; j < orig_bin_cnt; j++)
            sum += accum[i][j];
        
        double per_bin = sum / new_bin_cnt;

        unsigned int idx = 0;
        unsigned int bin = 0;
        
        while( idx < orig_bin_cnt )
        {
            double acc = 0;
            while( acc < per_bin && idx < orig_bin_cnt )
            {
                acc += accum[i][idx];
                new_bins[i][idx] = bin;
                idx++;
            }
            if( bin < new_bin_cnt - 1 )
                bin++;
        }
    }
}

void rebin_calls(const array_type5 &calls, const boost::multi_array<int, 2> &new_bins, unsigned int read_length, unsigned int orig_bin_cnt, unsigned int new_bin_cnt, array_type5 &new_calls)
{
    new_calls.resize(boost::extents[read_length][new_bin_cnt][4][4][4]);
    
    unsigned int i,j,ii,jj,kk;
    for(i = 0; i < read_length; i++)
        for(j = 0; j < orig_bin_cnt; j++)
            for(ii = 0; ii < 4; ii++)
                for(jj = 0; jj < 4; jj++)
                    for(kk = 0; kk < 4; kk++)
                        new_calls[i][new_bins[i][j]][ii][jj][kk] += calls[i][j][ii][jj][kk];
}

void convert_call_rates_to_probs(array_type5 &calls, unsigned int read_length, unsigned int bin_cnt)
{
    unsigned int i,j,ii,jj,kk;
    double sum;
    for(i = 0; i < read_length; i++)
        for(j = 0; j < bin_cnt; j++)
            for(ii = 0; ii < 4; ii++)
                for(jj = 0; jj < 4; jj++)
                {
                    sum = 0;
                    for(kk = 0; kk < 4; kk++)
                        sum += calls[i][j][ii][jj][kk];
                    
                    if(sum != 0)
                        for(kk = 0; kk < 4; kk++)
                            calls[i][j][ii][jj][kk] /= sum;
                }
}



/*********************************************************************************8
 * SEQUENCE Functions. 
*/


void rebin_seq_quals(const array_type4 &calls, unsigned int read_length, unsigned int orig_bin_cnt, unsigned int new_bin_cnt, boost::multi_array<int, 2> &new_bins)
{
    boost::multi_array<double, 2> accum(boost::extents[read_length][orig_bin_cnt]);
    new_bins.resize(boost::extents[read_length][orig_bin_cnt]);

    unsigned int i,j,ii,jj;

    for(i = 0; i < read_length; i++)
        for(j = 0; j < orig_bin_cnt; j++)
            for(ii = 0; ii < 4; ii++)
                for(jj = 0; jj < 4; jj++)
                    accum[i][j] += calls[i][j][ii][jj];
    
    for(i = 0; i < read_length; i++)
    {
        double sum = 0;
        for(j = 0; j < orig_bin_cnt; j++)
            sum += accum[i][j];
        
        double per_bin = sum / new_bin_cnt;

        unsigned int idx = 0;
        unsigned int bin = 0;
        
        while( idx < orig_bin_cnt )
        {
            double acc = 0;
            while( acc < per_bin && idx < orig_bin_cnt )
            {
                acc += accum[i][idx];
                new_bins[i][idx] = bin;
                idx++;
            }
            if( bin < new_bin_cnt - 1 )
                bin++;
        }
    }
}

// build a new call rates matrix from the old call rates and the new_bins.
void rebin_seq_calls(const array_type4 &calls, const boost::multi_array<int, 2> &new_bins, unsigned int read_length, unsigned int orig_bin_cnt, unsigned int new_bin_cnt, array_type4 &new_calls)
{
    new_calls.resize(boost::extents[read_length][new_bin_cnt][4][4]);
    
    unsigned int i,j,ii,jj;
    for(i = 0; i < read_length; i++)
        for(j = 0; j < orig_bin_cnt; j++)
            for(ii = 0; ii < 4; ii++)
                for(jj = 0; jj < 4; jj++)
                    new_calls[i][new_bins[i][j]][ii][jj] += calls[i][j][ii][jj];
}

// convert counts in call rates to probabilities.
void convert_seq_call_rates_to_probs(array_type4 &calls, unsigned int read_length, unsigned int bin_cnt)
{
    unsigned int i,j,ii,jj;
    double sum;
    for(i = 0; i < read_length; i++)
        for(j = 0; j < bin_cnt; j++)
            for(ii = 0; ii < 4; ii++)
            {
                sum = 0;
                for(jj = 0; jj < 4; jj++)
                    sum += calls[i][j][ii][jj];
                
                if(sum != 0)
                    for(jj = 0; jj < 4; jj++)
                        calls[i][j][ii][jj] /= sum;
            }
}
